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ABSTRACT 



Aims. The Gliese 581 planetary system has received attention because it has been proposed to host a low-mass planet in its habitable 
zone. We re-analyse the radial velocity measurements reported to contain six planetary signals to see whether these conclusions 
remain valid when the analyses are made using Bayesian tools instead of the common periodogram analyses. 

Methods. We analyse the combined radial velocity data set obtained using the HARPS and HIRES spectrographs using posterior 
sampling techniques and computation of the posterior probabilities of models with differing numbers of Keplerian signals. We do not 
fix the orbital eccentricities and stellar jitter to certain values but treat these as free parameters of our statistical models. Hence, we 
can take the uncertainties of these parameters into account when assessing the number of planetary signals present in the data, the 
point estimates of all of the model parameters, and the uncertainties of these parameters. 

Results. We conclude that based on the Bayesian model probabilities and the nature of the posterior densities of the different models, 
there is evidence in favour of four planets orbiting GJ 581. The HARPS and HIRES data do not imply the conclusion that there are 
two additional companions orbiting GJ 581. We also revise the orbital parameters of the four companions in the system. Especially, 
according to our results, the eccentricities of all the companions in the system are consistent with zero. 
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1. Introduction 

The nearby M dwarf Gliese 581 has received plenty of atten- 
tion during recent years. In 2005 a Neptune-mass planet can- 
didat e was found in its orbit using the HARPS spectrograph 
dBonfils et aU [2005). Two years later it was r eporte d to be a host 
to two additional super-Earths (lUdrv et all 120071) . In 2009 an 
new Earth-mass planet with a minimum mass of 1.9 M ffl was 
found in its orbit, making it a system with four pl anetary com- 
panions of relatively low mass dMavor et al.L 120091) . 

Ever since the discovery of the first companion orbiting GJ 
581, the star has been a target of intensive radial velocity (RV) 
surveys because few M dwarf stars are known to be hosts to plan- 
etary systems despite the fact that they are numerous even in the 
Solar neighbourhood. As a result, in 2010, GJ 581 was reported 
to have two more companions of planetary mass with GJ 581 
g, a 3 .1 M e planet, in the habitable zone of the star ( Vog t et ail 
[2010h . 

The discovery of this 6-planet system was made by analysing 
two high- precision RV data se ts made using the HARPS spec- 
trograph ilMavor et all 120091) and the HIRES spectrograph 
dVogt et all l2010h. The combined data set consists of 241 RV 
measurements. In lVogt et al.l(l2010h . the Keplerian signals of the 
six planets were discovered by studying the periodogram of the 
combined timeseries and by fitting the orbital parameters of the 
proposed companions to the data. 

The purpose of this letter is to see whether Bayesian 
data analysis gives consistent results with those reported by 



IVogtetalJ (120101) . Our major concern is, that by fixing eccen- 
tricities to zero the uncertainties of these eccentricities and their 
effect on the de tectability of plan etary signals were not taken 
into account by IVogt et alj (|2010). They did let the eccentrici- 
ties float freely and concluded that this did not produce any sig- 
nifican improvement to the fit. However, they did not discuss 
whether the uncertainty about the eccentricities could have an 
effect on the probabilities of finding the Kepler ian signals in the 
data in the first place. Also, in lVogt et al.l (120101) . the stellar jitter 
was estimated to have a value of 1.4 ms" 1 - simply because it 
yielded a reduced x 1 value of unity. Our second concern is that 
the uncertainty of the jitter was not taken into account either in 
the analyses. If its value was under- or overestimated, it could 
have a significant effect on th e detectability of the planetary sig - 
nals as demonstrated by e.g. iFordl (120061) : iGregorvl (l2007albl) : 
iTuomi & Kotirantal d2009l) . 
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In this letter, we reanalyse the combined RV data set of GJ 
581 using Bayesian tools - posterior samplings and model prob- 
abilities. First, we sample the parameter spaces of Keplerian 
models with k planetary companions by letting k - 0, ...,6. 
Second, we calculate the Bayesian probabilities for each of 
these models ik- Finally, we compare our results with those of 
IVogt et ahl (1201 Oh to see whether their periodogram analyses and 
fitting algorithms and Bayesian ones do yield similar results in 
this case. 
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2. Modelling and Bayesian model comparison 

2.1. Statistical model and posterior samplings 

Following Tuomi & Kotirantal d2009l) . we assume that the plan- 
ets do not interact with one another in the timescale of the mea- 
surements and model the superposition of k Keplerian signals 
simply by summing their effect on the RV. Consequently, there 
are five parameters describing the signature of an individual 
planet: RV amplitude (K), orbital period (P), orbital eccentric- 
ity (e), longitude of pericentre (oS), and the mean anomaly (Mo). 
As suggested by Ford (2006), to improve the efficiency of the 
sampling of the parameter space, we use the logarithm of the 
period in the samplings because it is a scale-invariant parameter. 

Our statistical model consists of the sum of Keplerian sig- 
nals and two sources of uncertainty. These sources are the instru- 
ment noise and the noise caused by the stellar surface - the stel- 
lar jitter. We model these as independent random variables with 
Gaussian density and zero mean. We assume that the standard 
deviation of the instrument no i se is known and use the values 
reported in lMavor etail (l2009l) : IVogt et all (l2010b . However, we 
adopt the standard deviation of the stellar jitter as a free param- 
eter of our statistical model and denote it as 07. The statistical 
model can be written as 



where the marginal integral is 



n.l = Zk(ti) +71 + €i + €j, 



(1) 



where ry is the measurement at time f,- made using telescope- 
instrument combination /, ik represents the k Keplerian signals, 
parameter j\ is the reference velocity, and e, and €j are Gaussian 
random variables describing the instrument noise, as reported 
by the observers, and the stellar jitter, respectively. As a result, 
there are 5 k + 3 free parameters in the parameter vectors 9k of 
our models - five parameters for each planet, jitter magnitude, 
and the reference velocities of the HARPS and HIRES measure- 
ments. 

We sample the parameter spaces of the different models us - 
ing the adaptive Metropolis algorithm of Haari o et"ai1 (12001 1) . 
This sampling algo rithm is similar to the famous Metrop olis - 
Hastings algorithm dMetropolis et al.Ul953l:lHast ings. 1973) but 
it adapts to the information gathered during the first n— 1 samples 
from the parameter posterior density by approximating this sam- 
ple as a multivariate Gaussian density. The nth sample is then 
drawn by using this multivariate Gaussian as a proposal density 
with the n - 1th value as a mean. While only an approximation, 
this algorithm converges to the posterior relatively rapidly - ap- 
parently even in the case of multimodal posterior sampled in this 
letter. 

When calculating the posterior densities for the parameters 
representing semimajor axes and RV masses of the planets, we 
took into account the uncertainty of the ste llar mass. This mass i s 
estimated to be 0.3 1 +0.02 M G for GJ 58 1 (iDelfosse et aU 12000)) . 
We sampled the densities of the semi-major axes and RV masses 
by drawing random numbers from the estimated density of the 
stellar mass that had a mean of 0.3 1 M and a standard deviation 
of 0.02 M Q . This enabled us to calculate more reliable estimates 
for the uncertainties of the semimajor axes and RV masses. 

2.2. Model comparisons 

We compare the models with differing number of planetary com- 
panions using the Bayesian model probabilities. The probability 
of the kth model is defined as 



P(r\z k )= f(r\ek,Zk)p(e k \zk)d0k, 



(3) 



P(zk\r) = 



P(r\z k )P(zk) 
Z p j=0 P(r\Zj)P(Zjy 



(2) 



f(r\6k,Zk) is the likelihood function, and p{9k\ik) the prior den- 
sity of the model parameters. In addition, P(ik) is the prior prob- 
ability of the A:th model and p is the greatest number of planetary 
signals in our analyses. 

We interpret the probabilities of Eq. (f2} as proper probabili- 
ties and require that the probability of confidently finding a kth 
planetary signal requires that P(ik) » P{ik-\)- In practice, ac- 
cording to Jeffreys (1961), we require that the probability of find- 
ing k signals is at least 150 times greater than that of finding k— 1 
signals to be able to claim confidently that there are k planets or- 
biting the target star. In addition, we require that the probability 
density of the kth planet has a unique maximum that can be in- 
terpreted as a Keplerian signal and is not likely to be caused by 
gaps in the data or by pure noise. 

Since the integral in Eq. (0 cannot be calculated analytically, 
we calculate its value from the sa mple from the posterior d ensity 
using the technique discussed in lChib & Jeliazkovl (1200 lh . 

3. Results 

We modelled the system by making two different assumptions 
about the orbital eccentricities of the planets. First, we treated the 
eccentricities as free parameters of the model, letting them vary 
freely in the sampli ngs of the param eter spaces. Second, follow- 
ing the analysis of Vogt et al.l d2010h . we tested a case where the 
eccentricities were only allowed to have low values - values be- 
low 0.2 - for the sake of dynamical stability of the system. 

3.1. Eccentricities as free parameters 

When adopting the orbital eccentricities of the k planets in the 
system as free parameters of the model, we receive the model 
probabilities (Pa) in Table Q] The probabilities of both model 
sets Pa and Pb, with different assumptions on the eccentricity, 
are scaled to unity according to Eq. (f2}. 

In TableQ] we present the probabilities of the different mod- 
els with k = 0, 6 planetary companions. We divide the model 
with k = 5 into two different solutions. In these solutions, 5a and 
5b are used to denot e the 5-companion models including the four 
planets reported bv lMavor et a l. (2009) and a fifth Keplerian sig- 
nal corresponding to one or other of the two signals proposed by 
IVogt et aLl d2010h with periods of roughly 37 and 433 days, re- 
spe ctively. The 6 pla net model includes all six planets reported 
by (jVogt et al., l2010h . Clearly, according to the probabilities, it 
cannot be concluded that there are the signals of six planetary 
companions in the data set. Instead, allowing the values of the 
eccentricities to be determined freely by the data leads to a con- 
clusion that there are clear signals of at least four companions 
in the data but the 4-companion model has such high probability 
that it cannot be ruled out confidently enough to claim that there 
are more than four Keplerian signals in the data. Hence, taking 
into account the uncertainties of the orbital and other parame- 
ters of the model s leads to results that contradict with those of 
IVogtet alJ<l2010h . 

The 5 -planet solution with eccentricities as free parameters 
consists broadly of two clear probability maxima for the 5th 
planet (the roughly 37 and 433 day periodicities) but we can- 
not conclude that the corresponding signals are real as opposed 
to artefact s of the da t a. Inte restingly, the 433 day periodicity pro- 
posed bv I Vogt et al.l d2010h has a greater probability than the 37 
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Table 1. Bayesian model probabilities of k planet models with 
free eccentricity (Pa) and eccentricity limited to values below 
0.2 (Pb)- The models 5a and 5b correspond to the solutions with 
the orbital period of GJ 581 f roughly at 37 and 433 days, re- 
spectively. 
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Fig. 1. The distribution of the orbital eccentricity of GJ 581 d 
from the 4-companion solution. A Gaussian curve with the same 
mean and variance is shown for comparison together with the 
mode, mean (//), standard deviation (<x), skewness (fj 3 ), and kur- 
tosis (jj 4 ) of the distribution. 



day periodicity but the former appears to consist of two closely 
spaced maxima in the periodicity space. There is one additional 
maxima in the vicinity of this period (at roughly 465 days), as 
can be seen in Fig. [2] 

If the orbital stability of the system is not considered, there is 
evidence in the data in favour of only fourplanetary companions. 
Curiously, we receive an interesting probability distribution for 
the eccentricity of planet GJ 581 d. Thi s distribution is sho wn in 
Fig.[TJand it supports the conclusion of May or etlril (12009 ), who 
claimed that this companion has an eccentricity of 0.38 ± 0.09. 
However, there is another maximum in this distribution close to 
zero, which makes the eccentricity of GJ 581 d consistent with 
zero. The eccentricities of the other companions were also con- 
sistent with zero, yet that o f GJ 581 c peaked a t 0.1 - a value 
consistent with the results of iMavor et al.l (120091) . 

According to our results, the point estimates of the orbital 
parameters and especially their uncertainty estimates need to be 
revised. The revised parameters and their 99% Bayesian credi - 
bility sets (2D0.99), as defined in e.g. lTuomi & Kotirantal J2009), 
are shown in Table [2] 

3.2. Low eccentricities 

We further assumed that the orbital eccentricities had values 
lower than 0.2 and repeated the analyses in the previous sub- 
section. The results of these analyses show that the five- and six- 
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Fig. 2. Same as Fig [TJ but for the period of GJ 581 f in the 5- 
companion solution. 



planet models do not have great enough posterior probabilities 
to be able to claim that there are more than fourplanetary com- 
panions orbiting GJ 581 (Table [TJ. T hese results show that the 
six-companion solution proposed by Vogtetal. (2010) cannot 
be considered to imply the existence of six planets in the system 
because the four-companion model cannot be shown to be an in- 
sufficient d escription of t he dat a confidently enough. Instead, the 
solution of Ma yor et"aD (120091) remains the most convincing so- 
lution even with the combined HARPS and HIRES dataset. Also, 
the updated parameter and uncertainty estimates of this solution 
are those presented in Table [2] 

Regarding the existence of the proposed companion in the 
habitable zone of GJ 581 with an orbital period of roughly 37 
days, the posterior probabilities in Table QJ imply that there is 
no evidence in favour of the existence of this companion. It is 
more probable that there is a companion corresponding to the 
periodicity of 433 days but even this periodicity appears to not 
be probable enough and also rather poorly constrained, as shown 
in Fig. [2] 

We also note, that limiting the eccentricities of the compan- 
ions to values lower t han 0.2 favours t he six-companion inter- 
pretation presented in Vogt e t al.l (l2010l) . The reason is that the 
eccentricity of GJ 581 d likely differs from zero and the result- 
ing solution where eccentricities are not limited describes the 
data much better than the limited case. 

Based on the data alone, the freely varying eccentricity is 
a much more likely scenario than the limited one with roughly 
25 times greater probability for the four-companion model and 
more than 100 times greater probability for the five-companion 
model. However, they are almost equal for the six-companion 
model, which means that the hypothetical six-companion model 
favours negligible eccentricities. 

3.3. Stellar jitter 

For the 4-companion solution, our estimate of stellar jitter of 
1.89 [1.51, 2.361 ms appears to differ significantly from the 
value reported bv lMavor et al.l (12009b of 1.2 ms -1 . First, we note 
that the jitter does not only correspond to the noise caused by 
the stellar surface, but it contains all the excess variations in 
the data not explained by the Keplerian model or the instrument 
noise. With measurements from two telescope-instrument com- 
binations, any small differences in the calibration of the tele- 
scopes and instruments may cause systematic differences to the 
measurements and lead to increased values for the jitter. Also, 
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Table 2. The four-planet solution of GJ 581 RV's. Maximum a posteriori estimates of the parameters and their ©0.99 sets. 



Parameter 


GJ 581 b 


GJ 581 c 


GJ 581 d 


GJ 581 e 


P [days] 


5.36849 [5.36810, 5.36888] 


12.916 [12.909, 12.922] 


66.85 [65.85, 67.76] 


3.1488 [3.1479, 3.1510] 


e 


0.007 [0, 0.051] 


0.09 [0, 0.29] 


0.39 [0, 0.67] 


0.08 [0, 0.43] 


A [Ills J 






1 m 8fi 9 971 

l.JO [U.OU, Z..Z.ZJ 


1 71 ri nfi 9 

L. / J [l.UU, Z..JJJ 


a> [rad] 


2.4 [0, 2tt] 


2.4 [0, 2tt] 


5.3 [0, 2tt] 


1.8 [0, 2tt] 


M [rad] 


1.4 [0, 2tt] 


3.1 [0, 2tt] 


4.9 [0, 2n\ 


1.6 [0, 2tt] 


m p sin i [M ffi ] 


15.70 [13.50, 17.70] 


5.53 [4.40, 6.93] 


4.51 [2.59, 6.69] 


1.81 [1.07, 2.55] 


a [AU] 


0.0405 [0.0380, 0.0430] 


0.0729 [0.0683, 0.0775] 


0.218 [0.203, 0.232] 


0.0283 [0.0267, 0.0303] 


y x [ms- 1 ] (HARPS) 


1.13 [0.48, 1.78] 








y 2 [ms^ 1 ] (HIRES) 


-0.35 [-0.95,0.19] 








a j [ms -1 ] 


1.89 [1.51, 2.36] 









undetected planets may i ncreas e the jitter between our results 
and those of lMavor et al.l (120091) because the observational time- 
line is longer in the combined dataset analysed here. 

Second, our lower limit of the jitter, or more accurately ex- 
cess noise, in the combined dataset is 1.51 ms" 1 - a value reason- 
ably close to the estimate of lMavor et al.l(l2009t) . The fact that we 
used posterior sampling technique, and took into account the un- 
certainties of all the parameters, can result in a greater estimate 
for the jitter magnitude because the common x 2 fitting, where 
jitter is not a free parameter, yields the lowest possible value for 
the jitter - fitting only the orbital parameters and reference ve- 
locities essentially corresponds to minimising the excess noise, 
not finding the most probable solution. 

3.4. Orbital stability 

We studied the orbital stability of out four-planet solution (Table 
|2]i briefly by using numerical integratio ns of the planetary orbits . 
We used the Bulirsch-Stoer algorithm (Bu lirsch & Stoerl [l966) 
because it is a relatively fast and reliable algorithm for studying 
the dynamics of planetary systems. To assess the instability of 
the system, we used the concept of Lagrange stability, in which 
the orbital parameters remain inside some bounded subspace of 
the parameter space for stable systems. Therefore, we considered 
the system not stable if there was orbital crossing, collision be- 
tween the planets, accretion of a planet by the star, or if a planet 
escaped from the system with a exceeding 100 AU. 

We drew a sample of 50 parameter vectors from the parame- 
ter posterior density by weighting the sample towards the high- 
eccentricity orbits - the orbital configurations most likely to be 
unstable - allowed by our solution. We integrated the orbits for 
10 000 years for each parameter vector because during that pe- 
riod of time even the outer planet would have completed more 
than 50 000 orbits. We checked whether the semi-major axes or 
orbital ecentricities evolved significantly during these integra- 
tions. 

In our the numerical integrations of the planetary orbits, the 
semi-major axes and eccentricities remained bounded in all but 
ten integrations. The semi-major axes remained almost constant 
and the orbital ecentricities of the three inner planets librated 
roughly between and 0.2, whereas the eccentricity of the outer 
planet librated only slightly around its initial value. Also, de- 
spite the possible moderate eccentricity of the outer planet of 
even more than 0.5, it did not appear to have a significant effect 
on the orbits of the three inner companions that would have re- 
sulted in orbital crossings. In ten cases, the initial eccentricities 
of planet c and e were greater than 0.2, where their probability 
densities have already become very low with respect to the MAP 
values. These eccentricities resulted in orbital crossings and it 



can be therefore concluded that these planets are likely to have 
eccentricities lower than 0.2. This is a consequence of the tight 
packing of the three innermost planets in the system. 

With these constraints, our sample from the posterior ap- 
peared to correspond to stable orbital configurations. Therefore, 
though having analysed the stability only briefly, we conclude 
that our revised solution likely corresponds to a stable system. 

4. Conclusions and discussion 

According to our analyses of the combined RV data for Gliese 
581 from HARPS and HIRES spectrographs, it cannot be con- 
cluded that there are six planetary companions orbiting Gliese 
581. We find that there are confidently four keplerian signals in 
the combined data set and that there is no evidence in favour 
of the existence of a low-mass planet in the habitable zone on 
Gliese 581 with an orbital period of 37 days. Therefore, we con- 
clude t hat our four-compan ion solution is consistent with the re- 
sults of lMavor et al.l d2009l) . though the uncertainties of the or- 
bital parameters and RV masse s require revision . We also con- 
clude that the interpretation of IVogt et al.l (1201 Ol) . who claimed 
that there are as many as six planets orbiting Gliese 581, appears 
to have been not supported by the data strongly enough. There is 
a periodicity of roughly 433 days in the combined data but it is 
not probable enough to claim that it is a real signal as opposed 
to an artefact of the data spacing or pure noise. 

We do realise that while the combined dataset provides evi- 
dence in favour of the existence of only four companions orbit- 
ing GJ 58 1 , it ca nnot be claimed that the proposed solution of 
Vogt etd] (120 101) is not real - it is fairly possible that the pro- 
posed companions f and g do exist in the system. However, the 
current data do not imply their existence in a statitically signifi- 
cant way, and the solution of the combined dataset remains that 
in Table [2] 

The most lik ely reasons for th e fact that we could not ver- 
ify the results of IVogt et al.l (1201 Ol) are that we i) treated the or- 
bital eccentricities of the planets as free parameters throughout 
the analyses; ii) adopted the magnitude of the stellar RV jitter 
as a free parameter of our statistical model instead of fixing its 
value to some a priori determined value or minimising it; and 
iii) calculated the posterior probabilities of the different models 
that take into account the Occamian principle of parsimony in a 
consistent manner - i.e. penalise the models more strongly the 
more free parameters they contain. We note that our estimate for 
the stellar j itter o f 1.89 ms~' contradicts the value obtained by 
IVogtetal.l(l2010j^ of 1.4 ms but is very close that estimated 
from the data of lWrighl (120051) of 1 .9 ms" 1 . 

According to the dynamical studies of our solution (Table|2), 
our revised solution is likely stable - a result consistent with the 
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analyses of iMavor et al.l (120091) . However, the availability of a 
sample from the posterior density of the orbital parameters al- 
lows a more detailed dynamical study of the Gliese 581 system. 
Such a study could help ruling out some subsets of the parame- 
ter space that appear to have a high probability based on the data 
alone and could not be shown to correspond to unstable config- 
urations in this study because of the relatively low integration 
times. 

Additional high-precision measurements of the low-mass 
target Gliese 581 are needed to assess the number of super- 
Earths in its orbit. Also, the first low-mass planet confidently or- 
biting its host-star within the limits of the stellar habitable zone 
remains to be confirmed with confidence. 
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